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Abstract — Proximity  ranging  is  very  important  in  many 
applications.  Ultrasonic  sensors  have  proven  to  be  very  cost 
effective  tools  for  this  purpose.  Most  of  the  currently  available 
time-of-flight  based  acoustic  proximity  ranging  systems  use 
the  conventional  matched  filter  based  time  delay  estimation 
approaches  to  measure  short  distances  between  proximity  objects. 
However,  in  the  presence  of  strong  and  closely  spaced  secondary 
echoes,  the  aforementioned  matched  filter  based  algorithms  tend 
to  fail  or  suffer  from  severe  performance  degradations  due  to  their 
poor  resolutions.  In  this  paper,  a  computationally  efficient  time 
delay  estimation  algorithm,  referred  to  as  multiecho  parameter 
estimation  for  acoustic  ranging  systems,  is  described  for  the  joint 
proximity  range  estimation  and  secondary  echo  mitigation.  The 
effectiveness  of  the  proposed  algorithm  is  demonstrated  by  both 
numerical  and  experimental  examples. 

Index  Terms — Acoustic  transducer,  proximity  ranging,  sec¬ 
ondary  echoes,  time  delay  estimation. 

I.  Introduction 

MEASURING  the  distance  between  a  known  base  loca¬ 
tion  and  the  surface  of  a  proximal  object  is  referred  to 
as  proximity  ranging  [1],  Proximity  ranging  is  very  important 
in  a  wide  range  of  remote  sensing  applications,  including  level 
detection,  robot  manipulation,  process  control,  nondestructive 
testing,  and  cavity  thickness  monitoring  [1] — [4], 

Various  types  of  sensors  based  on  different  physical  prin¬ 
ciples  such  as  capacitive  or  inductive  proximity  sensors,  laser 
displacement  sensors,  and  ultrasonic  sensors,  can  be  used  to 
perform  proximity  ranging  [3],  Among  these  sensors,  ultra¬ 
sonic  sensors  have  many  important  advantages  over  the  others. 
First,  they  can  operate  in  extreme  environmental  conditions, 
including  in  the  presence  of  fog,  dust,  dirt,  lighting  or  strong 
electromagnetic  interference  (EMI)  radiation.  Second,  ultra¬ 
sonic  sensors  can  be  used  for  accurate  distance  measurement 
by  using  the  time-of-flight  (TOF)  principle.  Moreover,  due  to 
the  low  speed  of  sound  in  air,  the  same  distance  resolution 
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can  be  achieved  by  using  much  simpler  electronic  circuits  for 
ultrasonic  sensors  than  those  for  laser  based  sensors.  Third, 
ultrasonic  sensors  can  be  fabricated  at  a  low  cost.  Since  gen¬ 
erating  ultrasound  only  requires  the  movement  of  a  surface 
and  can  be  effectively  implemented  by  using  electrostatic  or 
piezoelectric  devices,  the  fabrication  cost  for  ultrasonic  trans¬ 
ducers  is  usually  much  lower  than  that  of  laser  or  microwave 
based  sensors  [22],  Indeed,  ultrasonic  sensors  have  proven 
to  be  very  cost  effective  remote  sensing  tools  for  proximity 
ranging  [5]— [7], 

However,  due  to  the  presence  of  secondary  echoes,  using 
ultrasonic  sensors  for  very  close  proximity  ranging  is  very 
difficult.  When  ultrasonic  sensors  face  a  sound-hard  reflective 
surface,  the  reflected  sound  wave  can  bounce  back  and  forth 
several  times  between  the  sensors  and  the  reflection  surface 
before  decaying  to  zero,  which  results  in  unwanted  strong 
and  overlapping  secondary  echoes  in  the  received  signal.  The 
time  delays  for  the  secondary  echoes  are  approximately  integer 
multiples  of  the  time  delay  of  the  first  echo.  The  matched  filter 
based  methods  cannot  resolve  two  echoes  with  a  time  spacing 
less  than  the  reciprocal  of  the  signal  bandwidth  [8],  Hence, 
for  most  of  the  very  short  distance  measurement  scenarios, 
the  matched  filter  based  algorithms  tend  to  fail  or  suffer  from 
severe  performance  degradations  due  to  their  poor  resolutions 
[9],  Most  of  the  existing  super-resolution  time  delay  estimation 
approaches  (see  [10]  and  [11]  and  the  references  therein)  are 
developed  for  general  purposes.  They  do  not  exploit  the  a 
priori  information  of  the  integer  multiple  time  delays  and  the 
nonnegative  amplitude  due  to  the  acoustically  hard  reflections. 

In  this  paper,  we  focus  on  acoustic  proximity  ranging  in  the 
presence  of  secondary  echoes.  A  new  time  delay  estimation 
method  is  proposed  for  the  joint  proximity  ranging  and  sec¬ 
ondary  echo  mitigation.  We  establish  a  flexible  data  model  that 
takes  into  account  the  acoustically  hard  reflection  as  well  as 
small  random  time  delay  variations  of  the  secondary  echoes 
around  the  integer  multiples  of  the  time  delay  of  the  first 
echo.  The  new  time  delay  estimator  is  based  on  a  nonlinear 
least  squares  (NLS)  fitting  criterion.  However,  it  is  difficult  to 
directly  optimize  the  highly  nonlinear  NLS  cost  function.  We 
present  a  novel  computationally  efficient  time  delay  estimation 
algorithm,  referred  to  as  multiecho  parameter  estimation  for 
acoustic  ranging  systems  (multi-PEARS),  to  optimize  the  NLS 
cost  function. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  II,  we  describe  the  data  model  and  formulate  the 
problem  of  interest.  The  multi-PEARS  algorithm  is  presented 
in  Section  III.  Numerical  examples  and  experimental  results 


001 8-9456/03$  17.00  ©  2003  IEEE 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

APR  2003 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2003  to  00-00-2003 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Acoustic  Proximity  Ranging  in  the  Presence  of  Secondary  Echoes 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Florida, Department  of  Electrical  and  Computer 

Engineering, Gainesville, FL, 32611 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 

OF  PAGES 

13 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1594 


IEEE  TRANSACTIONS  ON  INSTRUMENTATION  AND  MEASUREMENT,  VOL.  52,  NO.  5,  OCTOBER  2003 


Fig.  1 .  Two-transducer  acoustic  proximity  ranging  system. 


are  given  in  Sections  IV  and  V,  respectively.  Section  VI 
concludes  this  paper.  The  derivations  of  the  single-echo  based 
PEARS  algorithm  (used  in  multi-PEARS)  and  the  Cramer-Rao 
bounds  (CRBs)  corresponding  to  multi-PEARS  are  included  in 
Appendices  A  and  B,  respectively. 

II.  Data  Model  and  Problem  Formulation 

As  shown  in  Fig.  1,  a  two-transducer  ranging  system  is  used 
to  measure  the  distance  between  the  transducers  and  an  acous¬ 
tically  hard  boundary.  We  assume  that  the  propagation  medium 
is  homogeneous  and  nondispersive  and  the  common  region  cov¬ 
ered  by  both  the  transmitter  and  the  receiver  is  a  smooth  planar 
surface.  The  transmitted  signal  is  emitted  from  the  transmitter 
toward  the  reflection  boundary.  The  directly  reflected  signal  re¬ 
ceived  by  the  receiver  is  called  the  first  echo.  Since  the  distance 
between  the  sensors  and  the  boundary  is  very  short,  the  sec¬ 
ondary  reflections  between  the  sensors  and  the  ranging  surface 
cannot  be  ignored.  The  time  delays  of  the  secondary  echoes  are 
approximately  integer  times  of  the  time  delay  of  the  first  echo. 
Hence  the  received  signal  can  be  modeled  as 
L 

y(t)  =  Y,  ais(t  —  It  -  Ai)  +  e(t),  0  <t<T  (1) 

i=i 

where  s(f),  0  <  /  <  T,  represents  the  known  (or  measured) 
real-valued  transmitted  signal  and  T  is  the  signal  duration;  y(t ) 
denotes  the  real-valued  received  signal;  [ 07 }  ,  denote  the  am¬ 

plitudes  of  the  received  signal  and  a\  is  assumed  to  be  nonneg¬ 
ative  due  to  the  acoustically  hard  reflection;  r  denotes  the  time 
delay  corresponding  to  the  first  echo,  and  the  time  delays  for 
the  secondary  echoes  are  It  +  Ai,  l  =  2, . . . ,  L,  with  Ai  =  0 
and  {A;}^_2  being  small  perturbations  around  the  integer  mul¬ 
tiples  of  r;  e(t )  is  the  real-valued  noise,  which  is  modeled  as  a 
zero-mean  Gaussian  random  process. 


The  received  signal  after  sampling  and  A/D  conversion  has 
the  form 

L 

y(nTs)  =  ^2  ais(nTs  -It-  A/)  +  e(nTs), 

1=1 

n  —  0.1,...:  A'  —  1  (2) 

where  Ts  denotes  the  sampling  period,  which  is  equal  to  the 
reciprocal  of  the  sampling  frequency  fs. 

Our  problem  of  interest  herein  is  to  estimate  r  from 
{y(nTs)}^Q  when  s(t),  0  <  t  <  T,  or  more  practically 
{s(n7Ta)}„=5)1  is  known. 

Instead  of  working  on  this  problem  in  the  time  domain,  the 
frequency  domain  is  preferred.  This  is  due  to  the  following 
two  reasons.  First,  for  the  time  domain  methods,  we  could 
either  restrict  r  to  be  multiples  of  Ts  or  resort  to  interpolation 
if  only  {s(nTg)}^L'j)1  is  known  [11].  This  inconvenience  can 
be  avoided  by  transforming  the  problem  into  the  frequency 
domain,  where  r  can  take  on  a  continuum  of  values.  Second, 
for  the  time  domain  data  model  in  (2),  {cq,A/}^L2  f°r  each 
secondary  echo  must  be  dealt  with  separately.  However,  by 
transforming  (2)  into  the  frequency  domain  and  using  the  facts 
that  {Ai}f_0  are  small  values,  a  more  concise  data  model  can 
be  obtained. 

Let  Y (to),  S(m),  and  E(m),  where  to.  =  —N/ 2,  —N/2  + 
1, . . . ,  N/2  —  1,  denote  the  discrete  Fourier  transforms  (DFTs) 
of  y(nTs),  s(nTs ),  and  e(nTs),  respectively.  Provided  that  the 
aliasing  due  to  sampling  is  negligible,  Y (to)  can  be  written  as 

L 

Y(m)  =  J2  <*iS(m)<N‘  (,r ' A,)m  +  E(m).  (3) 

1=1 

Since  both  the  transmitted  signal  s(t)  and  the  received  signal 
y(t)  are  real-valued,  their  Fourier  transforms  are  conjugate 


LI  et  al. :  ACOUSTIC  PROXIMITY  RANGING  IN  THE  PRESENCE  OF  SECONDARY  ECHOES 


1595 


symmetric .  Most  of  the  currently  available  ultrasonic  transducers 
are  resonant  devices  that  achieve  high  sensitivity  at  the  frequency 
of  less  than  200  kHz  and  are  narrowband  systems  [2].  Hence,  the 
transmitted  signal  s(t )  is  narrowband  in  general.  In  the  negative 
frequency  domain,  i.e.,  form  =  —N/2,  —N/ 2+1, . . . ,  0,  since 
s(t)  is  narrowband,  S(m )  mainly  occupies  a  few  frequency 
bins  in  the  range  [mo  —  Am,  rrio  +  Am],  where  mo  is  the 
peak  location  of  |S(m)|  and  2Am  ~  BSNTS  with  Bs  being 
the  signal  bandwidth.  Based  on  the  above  assumption  and  for 
very  small  { A; }  we  have  the  following  approximations 
for  the  signal  term  in  (3): 

aiS(m)<Cj^(lT+Al)m  = 

x  e~j^kAimo e-imal'rm 
M  aiS(m)e-jM  Arno^j^-lrm. 

=  frS(rn)e-j^lTm  (4) 


where  8i  =  aie~J^2^^NTa^Aim° ,  l  =  1,2 ,L.  Note  that 
8\  =  ft  |  is  nonnegative  and  {8/  }f=2  are  complex-valued.  In  the 
second  step  of  deriving  (4),  we  have  assumed  that  S(m)  ss  0 
for  m  ^  [mo  —  Am,  mo  +  Am]  and  {A/}£l2  are  so  small 
that  ((27r)/(iVTs))A/Am  ~  0  (i.e.,  ttAiBs  ~  0),  and  hence 
e-j(2n)/(NTs)Ai(m-m.0)  ^  l  for  m  £  [mo  —  Am, mo  +  Am]. 
Therefore,  the  small  time  delay  variations  of  the  secondary 
echoes  have  been  transformed  to  the  phases  of  their  complex 
amplitudes. 

Based  on  (3)  and  (4),  and  by  exploiting  the  conjugate  sym¬ 
metry  property  of  Y(m),  S(m),  and  E(m),  the  unknown 
parameters  can  be  obtained  by  minimizing  the  following  NLS 
criterion: 


A({/?,t})=  P2(m) 


Y(m)—S(m)^fiie  JjVTs 


1=1 


(5) 

where  f 3  =  [8\  82  •  •  •  8l]T  with  (-)T  denoting  the  transpose, 
{P(m)  =  1}~L-jv/2+i-  and  P(-N/ 2)  =  P( 0)  =  (l/v^)- 
Define 


x  = 


P(0)Y(0)]t 


Y 


where  K  =  N/2  +  1  is  the  dimension  of  the  vector  x.  Let  W 
be  a  K  X  K  diagonal  matrix 


W  =  diag{p(-f)S(-f). 

P  (“  2  _l)  5  (“  2  ■’)  •■  ' ' 

Let 


Define 


B(r)  =  [b(r)  b(2r)  •••  b(Pr)]  (10) 

and 

D(r)  =  WB(r).  (11) 

Then  (5)  can  be  expressed  as 

A({/?,t})  =  ||x-D(t)/?||2  (12) 

where  1 1  •  1 1  denotes  the  Euclidean  norm.  For  the  white  noise  case, 
the  above  NLS  approach  is  the  same  as  the  maximum  likelihood 
(ML)  method.  When  the  noise  is  colored,  the  NLS  approach  can 
still  give  very  accurate  parameter  estimates  [12], 

We  remark  that  the  time  delay  estimation  problem  considered 
herein  is  similar  to  the  well-known  harmonic  sinusoidal  param¬ 
eter  estimation  problem  [13] — [15].  However,  the  sinusoidal  sig¬ 
nals  in  our  problem  are  weighted  by  the  known  signal  spectrum 
and  the  amplitude  for  the  first  echo  is  nonnegative.  Since  most  of 
the  existing  harmonic  retrieval  methods  are  designed  for  signals 
with  complex-valued  amplitudes  and  flat  band-limited  spectra, 
they  are  not  directly  applicable  to  our  problem  of  interest. 

Note  that  minimizing  (12)  is  a  highly  nonlinear  optimization 
problem.  We  present  below  an  efficient  optimization  method, 
referred  to  as  the  multi-PEARS  algorithm,  to  minimize  (12)  by 
taking  into  account  both  nonnegative  81  and  the  weighted  har¬ 
monic  structure  of  the  data  model. 

III.  Multi-PEARS  Algorithm 

The  multi-PEARS  algorithm  optimizes  the  NLS  cost  func¬ 
tion  iteratively.  First,  we  assume  all  of  the  amplitudes  {81}^ 
are  complex-valued  and  obtain  an  initial  estimate  by  minimizing 
(12)  via  a  fast  Fourier  transform  (FFT)-based  method.  Next, 
based  on  the  initial  estimates,  we  subtract  out  the  components 
of  the  secondary  echoes  and  keep  only  the  component  of  the 
first  echo.  Then,  the  single-echo  based  PEARS  algorithm  [16] 
is  used  to  refine  the  initial  time  delay  estimate.  The  refining  step 
is  iterated  several  times. 

A.  Initialization  Stage 

Assume  that  {8i\i-i  are  complex-valued,  it  can  be  shown 
that  the  unknown  parameters  minimizing  (12)  can  be  deter¬ 
mined  by  [17] 


fo  =  argmaxco(r)  (13) 

T 

where 

c0(r)  =  xhD(t)  [Dh(t)D(t)]  1  D^(t )x  (14) 

and 


b(Zr) 


_ „•  2 -tv  It  ( N\  _  •  2  tv  It  ( W  1  1  \  "I 

e  - >  nts  \  2 )  e  J  NT»  1  2  1 )  ...  1 

[e-j2w/0Ir  e-j27r/i/r  _  .  .  e-j27r/K-i/rjT  (g) 


where 


fk  = 


1  k 

2TS  +  m3' 


k  =  0,1, . . . ,  K  —  1. 


(9) 


&=  [Dh(t)D(t)]  'D^rjxl^  (15) 

with  (-)H  denoting  the  conjugate  transpose.  According  to  the 
parsimony  principle  [18],  if  81  is  indeed  nonnegative,  the  es¬ 
timate  obtained  by  maximizing  co(r)  should  be  less  accurate 
than  that  obtained  by  minimizing  the  true  cost  function  of  (12). 
Note  that  (13)  is  a  one-dimensional  search  problem  and  cq(t) 
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usually  has  many  local  maxima.  Hence,  maximizing  co(r)  re¬ 
quires  the  search  over  a  very  fine  grid.  By  taking  advantage  of 
the  harmonic  structure  of  the  data  model,  we  use  zero-padding 
FFT  to  perform  such  a  search  computationally  efficiently. 

Note  that  DH(r)D(r)  in  f  14)  can  be  expressed  as 


D"(t)D(t)  = 


d0  d,\ 
d\ 


d 


L—l 


dL-1 


d1 

d\  d0 


(16) 


LxL 


to  the  nearest  integer.  Once  {di}j Iq1  and  {ai}[ J=1  are  in  hand, 
DH(r)D(r)  and  xHV(  t)  can  be  constructed  and  substituted 
into  (14)  to  obtain  the  cost  function  co(r)  at  the  r.  The  cal¬ 
culation  steps  for  the  initial  estimation  can  be  summarized  as 
follows. 

Step  1)  Zero-pad  {|FF(A,')|2}^q  and  {x* (k)W (&)}fc2!o  t0 
length  of  N  and  perform  N -point  FFTs.  Compen¬ 
sate  out  the  linear  phase  term  e-dN^)/ (2)  f0r  the  FFTs 
and  obtain  the  data  vectors  q  and  r,  respectively. 

Step  2)  For  a  given  r,  obtain  di  and  ai  as 


where  (•)*  denotes  the  complex  conjugate  and 


di  =  q 


K-l 


(. NtI ) 

(NT8) 


l  =  0, 1,  -  - ,  L 


(22) 


fc=0 


di  =  E  1^)1 

N 


2  e_-727r^,r 


=  eJ  2 

k= o 

Z  =0, 1, . . .  —  1 


2  e~^k 


and 

ai  =  r  ^ 

1 

_ i 

/  =  1, 2,  •  •  • ,  Z.  (23) 

,  , —  27 rZr 

NTS 

(17) 

Step  3)  Construct  Ds(r)D(r)  and  xffD(r)  by  using 
{di}fSQ  and  { a; obtained  in  Step  (2).  Compute 

with  W(k)  denoting  the  /,:th  diagonal  element  of  the  diagonal 
matrix  W  in  (7),  and  xHD(r)  can  be  expressed  as 


rH 


D(r)  =  [ai  a2 


ah] 


(18) 


where 


K-l 


ai=  [x*(k)W{k)]e~j2*fklT 
k= 0 

N_ 

=  ejl^J2  [x*(k)W(k)]  e~jujk 


k= 0 

l  =  1,2 


(19) 


Co(r)  for  the  given  r  based  on  (14). 

Step  4)  According  to  the  maximum  peak  location  of  the 
computed  Co(r)  in  Step  (3),  obtain  the  initial  esti¬ 
mate  to  of  t.  Once  fo  is  obtained,  0O  can  be  readily 
calculated  by  using  (15). 

In  practical  applications,  we  have  used  N  =  AN,  which  is 
sufficiently  accurate  for  the  initial  estimation  in  our  examples. 
Note  also  that  in  the  above  search  method,  we  need  to  compute 
the  FFTs  of  {|W/(fc)|2}^g  and  {x*  (k)W  (k)}^£ g  only  once. 
Thereafter,  the  look-up  table  method  is  used  to  find  {di  j/C",1 
and  { ai }fJ—1  for  a  given  r.  Thus,  this  method  is  computationally 
more  efficient  than  the  brute  force  search  method,  and  can  be 
easily  implemented  by  currently  available  dedicated  FFT  chips. 


with  x(k)  denoting  the  Atth  element  of  x  in  (6). 

It  can  be  seen  from  (17)  and  (19)  that  {d{\k~0  and  { ai }/i;1 
can  be  obtained  as  follows.  First,  the  discrete  time  Fourier 
transform  (DTFT),  which  can  be  efficiently  implemented 
by  using  zero-padding  FFT,  is  applied  to  the  sequences 
{|FF(A:)|2}^q  and  {x*(k)W(k) }^q,  respectively.  Next,  the 
linear  phase  term  g?'(tvW2)  js  multiplied  to  the  outputs  of 
DTFT,  and  finally,  the  sampling  is  performed  at  the  frequency 
points  lo  =  ( 2ttIt)/(NTs ). 

Define  the  following  two  vectors  which  are  obtained  by  mul¬ 
tiplying  the  phase  term  ej{Nu>)/(  2)  t0  the  FFTs  of{|FF(fc)|2}^2 
and  {x* (k)W (k)}^l_ g,  respectively,  as 


q=  k(o)  q(  1) 


q(N  —  1) 


(20) 


and 


r  =  r(0)  r(l) 


r(N  -  1) 


(21) 


where  N  denotes  the  data  length  after  zero-padding.  For  N  suf¬ 
ficient  large  and  for  any  given  r  and  /,  the  corresponding  di  and 
a/  can  be  approximately  obtained  by  picking  the  elements  of  q 
and  r  at  the  index  [(NIt)/(NTs)\,  where  [-J  denotes  rounding 


B.  Refining  Stage 

Assume  {0i}k_o  and  f  are  given,  we  subtract  out  the  sec¬ 
ondary  echoes  to  obtain 

L 

xi  =  x  -  ^^Wb(Zf).  (24) 

1-2 

Then  (12)  becomes 

C2  ({Pi, t})  =  ||xi  -/31Wb(r)||2  (25) 

where  is  assumed  to  be  nonnegative  due  to  the  acoustically 
hard  reflection.  We  have  developed  an  efficient  optimization 
algorithm,  referred  to  as  PEARS,  to  solve  the  time  delay  es¬ 
timation  problem  for  this  single  echo  case  [16].  To  make  this 
paper  self-contained,  the  derivation  of  the  PEARS  algorithm  is 
given  in  Appendix  A.  Note  that  PEARS  is  a  two-stage  algorithm 
which  first  obtains  an  initial  estimate  based  on  a  smooth  cost 
function,  and  then  refines  the  initial  estimate  based  on  the  true 
cost  function.  Herein,  since  we  have  already  obtained  the  initial 
estimate  of  r,  we  only  need  to  use  the  second  stage  of  PEARS  to 
refine  it.  Note  that  iteration  can  improve  the  accuracy  of  f .  The 
refined  estimate  f  can  be  used  to  redetermine  fl  via  (15).  Then 
the  refined  f  and  0  are  used  to  redetermine  xi  by  using  (24). 
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Fig.  2.  Numerical  results,  (a)  Transmitted  signal,  (b)  Received  signal  in 
the  absence  of  noise,  (c)  Comparison  of  the  RMSEs  of  the  estimates  and  the 
corresponding  CRBs  for  r.  Two  different  cases  for  the  time  delay  variations  of 
the  secondary  echoes  are  illustrated:  1)  A2  =  A3  =  0.  The  RMSEs  and  the 
corresponding  CRB  are  denoted  as  the  circles  and  the  solid  line,  respectively. 
2)  A 2  =  A 3  =  0.1  Te.  The  RMSEs  and  the  corresponding  CRB  are  denoted 
as  the  squares  and  the  dashed  line,  respectively. 


Based  on  the  updated  xi,  a  more  accurate  estimate  f  of  r  can 
be  obtained  by  using  PEARS.  With  the  above  simple  prepara¬ 
tions,  we  now  present  the  steps  of  the  multi-PEARS  algorithm. 
Step  1)  Obtain  the  initial  estimates  fo  and  0Q  by  using  the 
FFT -based  method  proposed  in  Section  III- A.  Let 
f  =  f0  and  0  =  0O. 

Step  2)  Obtain  xi  by  using  (24). 

Step  3)  Obtain  the  refined  f  by  using  the  second  stage  of 
PEARS  as  follows: 

Substep  3.1:  Obtain  the  measurement  for  each  fre¬ 
quency  point  as  xk  =  W*(k) x\(k), 

k  =  0, 1, -. . . ,  K  —  1,  where  x\ (k)  denotes 


the  A;th  element  of  xi . 

Obtain  the  phase  for  each  measurement  by 
Qk  =  —  arg(ffe),  k  =  0, 1, ...  ,7V  -  1, 
where  arg(x)  denotes  the  phase  of  x. 
Obtain  the  integers  zk  =  \_((2irfkT  — 
Ok)/ (27t))J ,k  =  0, 1, ... , K- 1,  where  fk 
is  defined  in  (9)  and  [-J  denotes  rounding 
to  the  nearest  integer. 

Obtain  the  refined  f  via 

t  =  (T,k=o  \xk\fk(Ok  +  2? TZk))/ 

Determine  refined  [1  via  (15)  by  using  the  refined  f . 
Iterate  Steps  2  through  4  and  calculate  the  residue 
after  each  iteration  as  ||x  —  D  (f)/3||2.  If  the  relative 
change  of  the  residues  between  two  consecutive  it¬ 
erations  is  less  than  a  small  value  e,  the  algorithm 
stops  (in  our  examples,  we  have  used  e  =  10-4  to 
test  the  convergence  of  the  algorithm). 

Before  we  proceed,  let  us  remark  on  the  following  two  facts. 
First,  multi-PEARS  converges  very  fast  due  to  the  good  initial 
estimates.  In  practice,  we  note  that  the  algorithm  usually  con¬ 
verges  within  three  or  four  iterations.  Second,  the  number  of 
the  echoes  L  is  determined  by  the  practical  measurement  envi¬ 
ronment.  Whenever  L  is  unknown,  it  can  be  estimated  from  x 
by  using,  for  instance,  the  generalized  Akaike  information  cri¬ 
terion  [12], 


Substep  3.2: 


Substep  3.3: 


Substep  3.4: 


Step  4: 
Step  5: 


IV.  Numerical  Examples 


To  demonstrate  the  estimation  accuracy  of  the  proposed 
multi-PEARS  algorithm,  we  present  a  numerical  example 
below  based  on  a  windowed  chirp  signal.  The  transmitted 
signal  is  defined  as 


s(t)  =  w(t )  cos 


27 vfct  +  7 


0  <  t  <  T 


(26) 

where  fc  denotes  the  carrier  frequency,  7  represents  the  chirp 
rate,  and 


wit)  = 


< 


0.5 -0.5  cos 

1, 

0.5 -0.5  cos 


TT(t-T0) 

T-u, 


0  <t<Tw, 

Tw  <  t  <  To  —  Tw, 

To  ~  Tw  <  t  <  To 

(27) 


with  Tw  =  To/10.  In  the  following  simulations,  we  use  N  = 
256,  7  =  2.57T  x  104,  the  signal  bandwidth  Bs  =  7 Tq/tt , 
the  carrier  frequency  fc  =  1 0/1,,  and  the  sampling  frequency 
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fs  =  32 Bs.  To  is  chosen  in  such  a  way  that  To  =  \/Ntt /647  = 
12.6  ms.  Thus  Ts  =  98.821  //,s,  Bs  =  0.316  kHz,  fc  = 
3.16  kHz,  fs  =  10.12  kHz,  and  T  =  25.3  ms.  Define  the  re¬ 
ciprocal  of  the  signal  bandwidth  as  Te=l/Bs  =  3.2  ms.  For 
the  received  signal,  r  =  16.2  Ts  =  0.506  Te,  ol\  =  0.5, 
q.2  =  0.3,  C13  =  0.1.  Note  that  the  time  delay  between  each 
echo  is  only  about  one  half  of  Te.  Thus,  the  first  echo  and  sec¬ 
ondary  echoes  cannot  to  be  resolved  by  using  the  matched  filter 
based  methods  in  this  case.  The  SNR  (signal-to-noise  ratio)  is 
defined  as  10  log w(a\Es/(Na2),  where  Es  is  the  energy  of  the 
signal  s(n,Ts).  The  root-mean-squared  errors  (RMSEs)  are  ob¬ 
tained  through  300  Monte  Carlo  trials. 

The  waveforms  of  the  transmitted  signal  and  the  noise  free 
received  signal  are  compared  in  Fig.  2(a)  and  (b).  In  Fig.  2(c), 
the  RMSEs  of  the  estimates  of  r  are  compared  with  the  cor¬ 
responding  CRBs,  which  are  the  best  performance  bounds  for 
any  unbiased  estimators  (the  CRBs  are  derived  in  Appendix  B). 
Two  different  cases  for  the  time  delay  variations  of  the  sec¬ 
ondary  echoes  are  illustrated:  1)  A2  =  A3  =  0,  and  2)  A2  = 
A3  =  0.1  Te.  In  case  (1),  the  RMSEs  of  the  estimates  and  the 
corresponding  CRB  are  denoted  as  the  circles  and  the  solid  line, 
respectively.  In  case  (2),  the  RMSEs  of  the  estimates  and  the 
corresponding  CRB  are  denoted  as  the  squares  and  the  dashed 
line,  respectively.  It  can  be  seen  that  when  there  are  no  time  delay 
variations  for  the  secondary  echoes  (i.e.,  A2  =  A3  =  0),  the 
RMSEs  of  the  estimates  can  achieve  the  corresponding  CRB. 
Note  also  the  threshold  effects  where  the  RMSEs  of  the  esti¬ 
mates  deviate  from  the  CRB  when  the  SNR  is  less  than  10  dB. 
When  there  exist  the  time  delay  variations  for  the  secondary 
echoes,  the  estimates  are  biased  and  the  RMSEs  of  the  estimates 
cannot  achieve  the  corresponding  CRB.  However,  if  the  time 
delay  variations  are  very  small  compared  with  the  reciprocal 
of  the  signal  bandwidth,  the  RMSEs  of  the  estimates  can  still 
approach  the  corresponding  CRB,  and  thus  our  method  can 
achieve  a  high  estimation  accuracy  in  these  cases  as  well.  In 
the  low  SNR  range,  the  RMSEs  of  the  estimates  deviate  from 
the  true  CRB  because  the  algorithm  can  no  longer  resolve  the 
true  peak  of  the  highly  oscillatory  cross  correlation  function 
between  the  received  signal  and  the  reference  signal.  However, 
in  this  case,  a  reliable  estimate  can  still  be  obtained  based  on 
the  peak  of  the  envelope  of  the  cross  correlation  function  [19], 
which  corresponds  to  obtaining  the  estimate  by  using  the  first 
stage  of  our  algorithm. 

Fig.  3(a)-(c)  illustrates  the  RMSEs  of  the  estimates  com¬ 
pared  with  the  CRBs  as  the  function  of  time  delay  by  fixing 
the  SNR.  The  time  delays  are  normalized  to  Te  and  we  assume 
A2  =  A3  =  0.  It  can  be  seen  that  as  the  time  delay  decreases, 
the  corresponding  CRB  increases.  This  is  as  expected  because 
the  smaller  the  time  delay,  the  more  the  echoes  are  overlapped, 
and  hence  the  worse  the  estimation  accuracy.  Noted  also  the 
threshold  effects  in  Fig.  3(a)  and  (b),  where  the  RMSEs  of  es¬ 
timates  deviate  from  the  CRB  when  the  time  delay  is  less  than 
a  certain  value.  In  order  to  still  achieve  the  corresponding  CRB 
in  the  smaller  time  delay  region,  the  higher  SNR  is  required. 

V.  Experimental  Results 

Experiments  have  been  conducted  by  using  the  commercially 
available  Panasonic  ceramic  transducers  (EFRRSB40K5  and 


Fig.  3.  RMSEs  and  CRBs  versus  the  normalized  time  delay  at  different  SNRs. 
The  time  delays  are  normalized  to  Te,  which  is  the  reciprocal  of  the  signal 
bandwidth,  (a)  SNR  =  10  dB,  (b)  SNR  =  15  dB,  and  (c)  SNR  =  20  dB. 
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TABLE  I 

Parameters  of  the  Panasonic  Ultrasonic  Ceramic  Transducers 


Part  No. 

EFRRSB40K5 

EFRTSB40K5 

Application 

Receiver 

Transmitter 

Resonant  Frequency  (kHz) 

40 

40 

Bandwidth  (kHz) 

4 

4 

Sensitivity  (dB)* 

-50 

— 

Sound  Pressure  Level  (dB)** 

— 

105 

Maximum  Input  Voltage  (Vrms) 

— 

20 

*  0  dB=l  V/Pa  **  0  dB— 2  x  IQ-5  Pa 


EFRTSB40K5).  Table  I  lists  the  parameters  of  the  transducers 
used  in  our  experiments  [20].  The  block  diagram  of  the  ex¬ 
perimental  setup  and  the  photograph  for  the  acoustic  ranging 
system  are  shown  in  Fig.  4(a)  and  (b),  respectively.  The  geo¬ 
metric  centers  of  the  two  transducers  are  17-mm  apart  and  they 
are  mounted  above  an  aluminum  plate  with  a  smooth  surface, 
which  is  used  to  simulate  the  planar  acoustically  hard  boundary. 
The  acoustic  transducers  are  fixed  on  a  three-dimensional 
(3-D)  traverse  controlled  by  a  personal  computer  (PC),  which 
is  used  to  adjust  the  position  of  the  transducers  related  to  the 
boundary  with  an  accuracy  of  0.4  //in.  The  excitation  signal 
for  the  transmitter  is  a  pulse  burst  generated  by  an  arbitrary 
waveform  generator.  The  pulse  repetition  frequency  is  192  Hz 
and  the  pulse  width  is  0.15  ms.  Each  pulse  has  been  modu¬ 
lated  by  a  carrier  frequency  of  40  kHz,  which  is  the  resonant 
frequency  of  the  transducer.  A  16-bit  ADC  is  used  to  sample 
( fs  =  196608  Hz)  the  excitation  signal  for  the  transmitter 
and  the  received  signal  from  the  receiver.  The  recorded  data 
are  transferred  to  a  host  PC  and  the  multi-PEARS  algorithm  is 
applied  to  the  recorded  data  to  determine  the  distance  estimates. 
The  sound  speed  is  obtained  by  [21] 

c  =  sJ-fBTr  (28) 

where  7  =  1.4  is  the  ratio  of  specific  heats,  R  =  287  J / (kg  •  K) 
is  the  gas  constant,  and  Tr  =  (273.16  +  21)  K  in  our  experi¬ 
ments.  Thus  we  have  c  =  344  m/s.  In  our  experiments,  since 
the  transmitter  and  the  receiver  are  closely  placed,  there  exists  a 
small  amount  of  crosstalk.  Note  that  the  crosstalk  is  unchanged 
and  independent  of  the  reflected  signal,  and  hence  it  can  be 
recorded  and  stored  in  the  system  in  advance  and  subtracted 
out  from  each  received  signal  thereafter.  All  the  received  sig¬ 
nals  below  are  illustrated  after  the  crosstalk  subtraction. 

Fig.  5  shows  the  recorded  signals  in  the  experiments. 
Fig.  5(a)  is  the  transmitted  signal,  which  is  obtained  when 
the  transmitter  and  the  receiver  are  face  to  face.  The  distance 
between  the  transducers  has  been  measured  precisely.  Thus 
this  waveform  is  taken  as  the  known  transmitted  waveform. 
Fig.  5(b)  shows  the  discrete  Fourier  transform  (magnitude) 
of  the  signal  in  Fig.  5(a).  It  can  be  seen  that  the  transmitted 
signal  is  narrowband  with  the  bandwidth  about  4  kHz.  Thus, 
the  nominal  minimum  distance  required  to  resolve  the  first  and 
secondary  echoes  by  the  matched  filter  method  is  about  43  mm. 

We  then  increase  the  distance  between  the  transducers  to  the 
boundary  starting  at  29  mm.  The  entire  measurement  range  is 
exactly  10  mm  and  is  covered  in  21  discrete  levels.  The  distance 
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Fig.  4.  Experimental  setup  for  the  acoustic  ranging  system,  (a)  Block  diagram 
of  the  experimental  setup,  (b)  Photograph  of  the  setup  for  the  sensors  and  the 
boundary. 


difference  between  two  consecutive  levels  is  0.5  mm.  Fig. 5(c) 
to  (f)  show  the  measured  waveforms  recorded  at  the  levels  3,  5, 
15,  and  19,  respectively.  As  we  can  see  that  due  to  the  presence 
of  the  secondary  echoes,  the  shapes  of  the  received  signals  at 
different  levels  are  quite  different.  For  the  noise  in  this  experi¬ 
ment,  the  electronic  noise  dominates  and  the  SNR  is  about  20  dB 
based  on  the  definition  in  Section  IV.  In  the  practical  environ¬ 
ment,  the  noise  level  can  be  much  higher  due  to  the  presence  of 
the  environmental  noise.  However,  for  a  steady  and  controllable 
environmental  noise,  one  can  obtain  it  by  using  another  receiver 
independently  driven,  and  remove  it  from  the  real  received  sig¬ 
nals  thereafter. 

Fig.  6  shows  an  example  where  the  multi-PEARS  algorithm 
is  applied  to  the  recorded  signal  at  level  19.  We  use  generalized 
Akaike  information  criterion  to  determine  the  number  of  the 
echoes  and  obtain  L  =  4  for  this  example.  It  can  be  seen  from 
this  figure  that  multi-PEARS  can  correctly  resolve  the  first  and 
secondary  echoes  even  though  they  are  heavily  overlapped.  Note 
also  that  the  reconstructed  signal  matches  the  original  signal 
very  well.  The  relative  error  between  these  two  signals  is  only 
0.7%. 

Fig.  7  shows  the  measurement  results  for  the  entire  range 
with  the  different  choices  for  L  (the  number  of  echoes).  If 
the  secondary  echoes  are  not  taken  into  account  and  L  =  1 
is  used,  an  accurate  f  cannot  be  obtained.  By  using  L  =  2, 
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Fig.  5.  Recored  waveforms  in  the  real  experiment,  (a)  Transmitted  waveform,  (b)  Discrete  Fourier  transform  of  signal  in  (a),  (c)  Received  signal  at  level  3. 
(d)  Received  signal  at  level  5.  (e)  Received  signal  at  level  15.  (f)  Received  signal  at  level  19. 


the  measurement  results  are  improved.  But  there  are  still  four 
levels  where  the  measurements  are  quite  different  from  the 


true  distances.  (The  results  for  L  =  3  are  similar  to  those  of 
L  =  2.)  By  using  L  —  4,  very  good  results  are  obtained  for 
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Fig.  6.  Application  of  multi-PEARS  to  a  recorded  signal  at  level  19.  (a)  The  original  signal,  (b)-(e)  The  reconstructed  signals  by  using  multi-PEARS  estimates 
for  the  echoes  1  to  4,  respectively,  (d)  The  reconstructed  signal  by  superposition  of  the  signals  in  (b)-(e). 


all  distances.  The  zoomed-in  view  of  the  measurement  errors  error  is  about  0.02  mm  for  L  —  4.  The  ratio  between  this 
are  shown  in  Fig.  7(b).  Note  that  the  maximum  measurement  error  and  the  entire  measurement  range  is  2  X  10-3. 
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Fig.  7.  Stationary  level  measurement  results  obtained  by  applying  the 
multi-PEARS  algorithm  to  the  measured  signal  with  the  different  choices  of 
L.  The  ranging  distance  covers  from  29  to  39  mm.  The  distance  between  two 
consecutive  levels  is  0.5  mm.  (a)  Measurement  results  versus  the  true  distance 
obtained  by  traverse,  (b)  Measurement  errors. 


VI.  Conclusion 

Due  to  the  presence  of  unwanted  strong  and  closely  spaced 
secondary  echoes,  it  is  very  difficult  to  obtain  accurate  short 
distance  measurements  for  acoustic  proximity  ranging  systems 
using  traditional  matched  filter  based  methods.  In  this  paper, 
a  computationally  efficient  time  delay  estimation  algorithm, 
referred  to  as  multi-PEARS,  is  presented  for  the  joint  proximity 
ranging  and  secondary  echo  mitigation.  Numerical  results 
demonstrate  that  multi-PEARS  can  deal  with  secondary  echoes 
effectively,  and  the  estimates  obtained  by  using  multi-PEARS 
can  approach  the  corresponding  CRBs  as  the  SNR  increases. 
Experimental  results  obtained  by  using  Panasonic  ceramic 
transducers  show  that  multi-PEARS  can  also  perform  very  well 
in  a  practical  environment. 


Appendix  A 
Derivation  of  PEARS 

Based  on  the  same  notation  for  (25),  the  NLS  cost  function 
for  the  single  echo  case  is 

Co({(3i ,r})  =  ||Xl  -/31Wb(r)||2  (29) 

where  ||  •  ||  denotes  the  Euclidean  norm.  Note  that  is  nonneg¬ 
ative.  Hence  the  cost  function  £2  in  (29)  can  be  expressed  as 

£2  =  [xi  -  /3iWb(r)]ff  [xi  -  ftWb(r)] 

=  xfxi  -  ^ixfWb(r)  -  /Jib^fVjW^xi 
+  ft2bff(r)WJWb(r) 

=  xfxi  +  f3fbH(T)WHWb(r) 

-  2/3iRe  [bH(r)WHxi]  (30) 

where  Re '(;/;)  is  the  real  part  of  x.  Since  xf*xi  + 
fibH(r)WHWb  (t)  is  a  constant  and  j3i  is  nonnega¬ 
tive,  the  r  that  minimizes  £2  is  the  one  that  maximizes  the  term 
2/3iRe[bH(r)WHxi],  which  has  the  form 

fc  =  argmaxci(r)  =  argmaxRe  [b^MW^xi]  .  (31) 

Note  that  the  cost  function  ci(r)  can  be  expressed  as 

Ci(r)  =Re  [bH(r)WHxi] 

K-i 

=  ^2  l^fcl  cos(2-7r/feT  —  Ok)  (32) 

fc= 0 


where 


Xk  =  W*(k)xi(k),  k  =  0,l,...,K  -1  (33) 

and 

Ok  =  ~  arg (xk),  k  =  0, 1, . . . ,  K  -  1  (34) 

with  arg(a;)  being  the  phase  of  x.  In  the  absence  of  noise, 

8k  =  2irfkT  -  2nzk,  k  =  0, 1,  1  (35) 

with  {zk}kl 2  t>eing  some  integers.  It  can  be  seen  from  (32)  that 
ci(t)  is  a  sum  of  K  cosine  functions  of  r  with  frequencies  /q 
to  fh'-  \  >  and  hence  is  highly  oscillatory,  which  makes  it  very 
difficult  to  find  the  global  maximum  directly. 

If  l3]  were  complex-valued,  it  can  be  shown  that  minimizing 
£ 2  with  respect  to  r  would  yield 

fC2  =  argmaxc2(r)  =  argmax lbH(r)WHxi| .  (36) 

Since  C2(r)  is  much  smoother  than  the  true  cost  function  Ci(r), 
we  can  obtain  an  initial  estimate  based  on  C2(t)  and  then  re¬ 
fine  it  based  on  c,\  (r) .  The  first  stage  of  PEARS  is  to  deter¬ 
mine  an  initial  estimate  based  on  C2(r).  Note  that  02(7")  can 
be  efficiently  maximized  by  applying  FFT  with  zero-padding 
to  Wflxi.  Once  the  initial  estimate,  say  fo,  is  available,  we  can 
refine  it  based  on  the  true  cost  function  c-\  (r) . 
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Differentiating  Ci(r)  with  respect  to  r  yields 
K—l 

c'i(t )  =  -^2  l^fclWfc  sin(27r/fcT  -  6k).  (37) 

k= 0 

The  goal  here  is  to  approximate  (37)  around  fo  with  the  first- 
order  approximation  to  find  an  approximate  solution  to  <\  (r)  = 
0,  which  corresponds  to  the  local  maximum  of  Ci(r)  around  fo- 
From  (35),  we  have  the  following  approximations  at  the  high 
SNR: 

27r/fcfo  ss  2irzk  +  8k,  k  =  0, 1, . . . ,  K  -  1  (38) 


zk=  2irfkl°  °k  ,  k  =  0,1, . . .  ,K  —  1  (39) 

Z7T 

with  L-J  denoting  rounding  to  the  nearest  integer.  In  a  small 
range  around  fo,  using  (38)  yields  the  following  approxima- 


sin(27r/fcT  -  8k)  =  sin(27r fkT  -  6k  -  2tt zk) 

~  sm[27r/fc(r-f0)] 

~  2nfkT  -  6k  -  2nzk.  (40) 

Therefore,  c^(r)  around  fo  can  be  approximately  expressed  as 
K- 1 

4 (r)  ~  -  ^2  27r|a;fc|/fc(27r/fcr  -  I %  -  2irzk).  (41) 
k= 0 

Simplifying  (41)  and  setting  it  to  zero,  we  obtain  the  refined 
solution  of  r  as 

|  =  EfrJo1 1 Xk\hih  +  2? xzk) 

>:k~:  24,',!./- 

This  estimator  can  be  interpreted  as  a  weighted  sum  of  K  esti¬ 
mates  of  r  obtained  from  frequencies  /q  to  Jk-\  separately. 


Er  =  ]E  2-J  E( 0) 

s"=K4+i)  5H+2 

•••  S(-l)f 


sr  =  |5  5(0)  j  (48) 

g(/r  -  A/)  —  e~20fVT+Ak>(~% ) 

f;  ■' (,r  •  A,)(  f  •  |)  ...  1  T(49) 

gc(/r  +  A,)  =  e-jT$k{lT+Al)(-%+1} 

0t+ A; )  (-  4-  +2) 

e~i T  (50) 

gr(lr  +  Ai)  =  )  1  (51) 

G  =  g(r)  g(2r  +  A2) 

•  •  •  g  (Lt  +  Ai)]  (52) 

Gc  =  [gc(f)  gc(2r  +  A2) 

•  •  •  g  c(Lt  +  Al)]  (53) 


Gr  =  [gr(r)  g,.(2r  +  A2) 


g r(Lr  +  Ai)] .  (54) 


We  have 


Yc  —  nc  +  Ec 


Yr  —  Mr  4-  E  r 


Me  —  ScGco: 


Appendix  B 


Derivation  of  the  CRBs 

We  sketch  below  the  derivation  of  the  CRBs  for  the 
data  model  in  (3)  where  {a.]  }\J=\  are  assumed  to  be 
real-valued.  Due  to  the  conjugate  symmetry  property  of 
DFTs,  {Y {my\2,j=-(N / 2)  can  exPressed  in  terms  of 
{^(™)}™=_(7v/2)  with  {Y(m)}“1=_(iv/2)+1  being  com¬ 
plex-valued  and  {Y (—(N/2)),Y (0)}  being  real-valued.  Let 


E(-  1)]T  (45) 


Mr  =  S,.Gra  (58) 

with 

ot  =  [ot  1  a2  ■■■  aL]T.  (59) 

Assume  that  the  additive  noise  {e(?tTs)}^_r01  is  a  real-valued 
zero-mean  white  Gaussian  random  process  with  variance  a2. 
Denote 

6  =  [a  1  •••  aL  t  A2  •••  Ai  cr2]T .  (60) 

Since  the  DFT  matrix  is  a  unitary  operator,  the  joint  probability 
density  function  for  {Yc,  Y,,}  has  the  form 

p({Yc,Yr}|0)=  \  expj-^(Yc-Mc)g(Yc-Mc) 

(27t)"2"  aN  {  <J~ 

~  7^ (Yr -Mr)T(Yr -Mr) |  •  (61) 
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Based  on  the  extended  Slepian-Bang’s  formula,  the  CRB  ma¬ 
trix  can  be  computed  as 


[CRB-1(0)] ..  =  -tr  [Q-LQ^Q-'Q,.;.]  +  [^.f  Qr“V,' 


xtr  [Qc  lQCj]  +  2  Re 


rj 

(62) 


where 

dg(/r  +  AQ 
dAi 


=  ~J 


2ir 

NT. " 


where  Q,,  =  <j2Ir  with  I,,  being  the  2x2  identity  matrix, 
Qc  =  cr2Ic  with  Ic  being  the  (N / 2  —  1)  X  (N / 2  —  1)  identity 
matrix,  and  X'  denotes  the  derivative  of  X  with  respect  to  the 
/th  unknown  parameter.  Since  Qr  and  Qc  do  not  depend  on  the 
parameters  in  //,,  and  fj,c,  and  ftr  and  (ir  do  not  depend  on  the 
elements  in  Q,,  and  Qc,  it  can  be  shown  that  the  matrix  CRB(0) 
is  block  diagonal  with  its  last  row  and  last  column  being  zero 
except  for  the  last  diagonal  element.  Let  the  signal  parameter 
vector  rj  be  denoted  as 


and 


F  = 


Then 


dn 

da\ 


J  NTc 


1  T 


(/t+AL(—  iy+l) 


dn 

dar , 


On 

8t 


dn 

dA o 


CRB(r?)  =  [Re(FffF)]‘ 
which  is  easily  evaluated  numerically. 


0 


dn 


(71) 


(72) 


(73) 


V  =  [cm  •  •  •  aL  t  A2  •  •  •  Al]t.  (63) 


References 


Let 


H  =  WGa 


(64) 


with  W  being  defined  in  (7),  and  Q  =  cr2l  with  I  being  the 
identity  matrix  of  dimension  iV/2  +  1.  Then 


[CRB-1  (77)]^.  =  M,fQr-V,' 


+  2  Re 


Mci  Qc  M, 


1..  / 

Cf 


=  2  Re 


,4hQ  V* 


(65) 


(68) 


with 

dg{h  +  A[) 
8t 


=  ~J 


.  2irl 


NT ' 


or 


(69) 


[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 


Define 

[9] 

g=wg(/T  +  i,), 

l  =  1, 2  (66) 

[10] 

dr  Or 

(67) 

[ii] 

where 

[12] 

<9G 

'Sg(r)  3g(2  r  +  A2) 

8g(Lr  +  AL)' 

[13] 

7h  ~ 

8t  8t 

8t 

[14] 

[15] 

[16] 

[17] 

[18] 
[19] 


and 


d (i 

dAi 


=  W 


dg  (It  +  A/) 
dA  1 


at,  /  =  2, 3, . . .  ,1/  (70) 


[20] 

[21] 

[22] 


W.  Gopel,  J.  Hesse,  and  J.  N.  Zemel,  Sensors:  A  Comprehensive 
Survey.  Weinheim,  Germany:  Wiley- VCH,  1989. 

O.  Brand  and  H.  Baltes,  “Micromachined  resonant  sensors:  an 
overview,”  in  Sensors  Update,  H.  Baltes,  W.  Gopel,  and  J.  Hesse, 
Eds.  Weinheim,  Germany:  Wiley-VCH,  1998,  vol.  4,  pp.  3-51. 

W.  Manthey,  N.  Kroemer,  and  V.  Magori,  “Ultrasonic  transducers  and 
transducer  arrays  for  applications  in  air,”  Meas.  Sci.  Technol.,  vol.  3,  pp. 
249-261,  1992. 

S.  Ashley,  “Wrap  drive  underwater,”  Sci.  Amer.,  vol.  284,  pp.  70-79, 
May  2001. 

D.  P.  Massa,  “Choosing  an  ultrasonic  sensor  for  proximity  or  distance 
measurement.  Part  1,”  Sensors,  vol.  16,  pp.  34-37,  Feb.  1999. 

Massa  Products  Corporation.  [Online]  http://www.massa.com. 
Migatron,  Inc.  [Online]  http://www.migatron.com. 

R.  Wu  and  J.  Li,  “An  efficient  algorithm  for  time  delay  estimation,” 
IEEE  Trans.  Signal  Processing,  vol.  46,  pp.  2231-2235,  Aug.  1998. 

J.  E.  Ehrenberg,  T.  E.  Ewart,  and  R.  D.  Morris,  “Signal-processing  tech¬ 
niques  for  resolving  individual  pulses  in  a  multipath  signal,”  J.  Acoust. 
Soc.  Amer.,  vol.  63,  pp.  1861-1865,  June  1978. 

R.  Wu  and  J.  Li,  “Time-delay  estimation  via  optimizing  highly  oscilla¬ 
tory  cost  functions,”  IEEE  J.  Ocean.  Eng.,  vol.  23,  pp.  235-244,  July 
1998. 

R.  Wu,  J.  Li,  and  Z.  Liu,  “Super  resolution  time  delay  estimation  via 
MODE-WRELAX,”  IEEE  Trans.  Aerosp.  Electron.  Syst.,  vol.  35,  pp. 
294-307,  Jan.  1999. 

J.  Li  and  P.  Stoica,  “Efficient  mixed- spectrum  estimation  with  applica¬ 
tions  to  target  feature  extraction,”  IEEE  Trans.  Signal  Processing,  vol. 
44,  pp.  281-295,  Feb.  1996. 

H.  Li,  P.  Stoica,  and  J.  Li,  “Computationally  efficient  parameter  esti¬ 
mation  for  harmonic  sinusoidal  signals,”  Signal  Processing,  vol.  80,  pp. 
1937-1944,  2000. 

B.  G.  Quinn  and  P.  J.  Thomson,  “Estimating  the  frequency  of  a  periodic 
function,”  Biometrika,  vol.  78,  pp.  65-74,  Mar.  1991. 

M.  Zeytinoglu  and  K.  M.  Wong,  “Detection  of  harmonic  sets,”  IEEE 
Trans.  Signal  Processing,  vol.  43,  pp.  2618-2630,  Nov.  1995. 

X.  Li,  R.  Wu,  S.  Rasmi,  J.  Li,  L.  N.  Cattafesta,  and  M.  Sheplak,  “An 
acoustic  proximity  ranging  system  for  monitoring  the  cavity  thickness,” 
IEEE  Trans.  Ultrason.,  Ferroelect.,  Freq.  Contr.,  vol.  50,  pp.  898-910, 
July  2003. 

P.  Stoica  and  R.  L.  Moses,  Introduction  to  Spectral  Analysis.  Upper 
Saddle  River,  NJ:  Prentice-Hall,  1997. 

T.  Soderstrom  and  P.  Stoica,  System  Identification.  London,  U.K.: 
Prentice-Hall  International,  1989. 

X.  Li,  R.  Wu,  M.  Sheplak,  and  J.  Li,  “Multifrequency  CW-based  time- 
delay  estimation  for  proximity  ultrasonic  sensors,”  Proc.  IEEE  Radar, 
Sonar  and  Navigation,  vol.  149,  pp.  53-59,  Apr.  2002. 

Panasonic  Corp.,  Ultrasonic  Ceramic  Sensors,  Burlington,  MA,  2001. 

J.  D.  Anderson,  Modern  Compressible  Flow.  New  York:  McGraw- 
Hill,  1990. 

J.  Jacob  Fraden,  Handbook  of  Modern  Sensors,  2nd  ed.  Woodbury, 
New  York:  Amer.  Inst.  Phys.,  1997,  ch.  5,  pp.  233-278. 


LI  et  ah'.  ACOUSTIC  PROXIMITY  RANGING  IN  THE  PRESENCE  OF  SECONDARY  ECHOES 


1605 


Xi  Li  (S’ 01)  received  the  B.Sc  and  Ph.D.  degrees  in 
electrical  engineering  from  Nanjing  University  of 
Science  and  Technology  (NUST),  Nanjing,  China,  in 
1995  and  1999,  respectively.  He  is  currently  pursuing 
the  Ph.D.  degree  in  electrical  engineering,  with  a 
minor  in  aerospace  engineering,  at  the  University  of 
Florida,  Gainesville. 

Since  May  2000,  he  has  been  a  Research  Assistant 
in  the  Department  of  Electrical  and  Computer  En¬ 
gineering  at  the  University  of  Florida.  His  research 
interests  include  spectral  estimation  and  signal  pro¬ 
cessing  for  acoustic  and  radar  applications. 


Jian  Li  (S’87-M’91-SM’97)  received  the  M.Sc. 
and  Ph.D.  degrees  in  electrical  engineering  from  The 
Ohio  State  University,  Columbus,  in  1987  and  1991, 
respectively. 

From  April  1991  to  June  1991,  she  was  an  Adjunct 
Assistant  Professor  with  the  Department  of  Electrical 
Engineering,  The  Ohio  State  University.  From  July 
1991  to  June  1993,  she  was  an  Assistant  Professor 
with  the  Department  of  Electrical  Engineering,  Uni¬ 
versity  of  Kentucky,  Lexington.  Since  August  1993, 
she  has  been  with  the  Department  of  Electrical  and 
Computer  Engineering,  University  of  Florida,  Gainesville,  where  she  is  cur¬ 
rently  a  Professor.  Her  current  research  interests  include  spectral  estimation, 
array  signal  processing,  and  signal  processing  for  wireless  communications  and 
radar.  She  is  currently  a  Guest  Editor  for  Multidimensional  Systems  and  Signal 
Processing. 

Dr.  Li  is  a  member  of  Sigma  Xi  and  Phi  Kappa  Phi.  She  received  the  1994 
National  Science  Foundation  Young  Investigator  Award  and  the  1996  Office  of 
Naval  Research  Young  Investigator  Award.  She  is  currently  an  Associate  Editor 
for  the  IEEE  Transactions  on  Signal  Processing. 


Renbiao  Wu  (M’95-SM’01)  received  the  B.Sc.  and 
M.Sc.  degrees  in  electrical  engineering  from  North¬ 
western  Polytechnic  University,  Xian,  China,  in  1988 
and  1991,  respectively,  and  the  Ph.D.  degree  in  elec¬ 
trical  engineering  from  Xidian  University,  Xian,  in 
1994. 

From  May  1994  to  February  1996,  he  was  a 
Postdoctoral  Fellow  at  the  College  of  Marine 
Engineering,  Northwestern  Polytechnic  University, 
where  he  was  promoted  to  Associate  Professor  in 
December  1995.  From  March  1996  to  February 
1997,  he  was  a  Visiting  Scholar  at  the  Center  for  Transportation  Research, 
Virginia  Polytechnic  Institute  and  State  University,  Blacksburg.  From  March 
1997  to  December  1998,  he  was  a  Visiting  Scholar  in  the  Department  of 
Electrical  and  Computer  Engineering,  University  of  Florida,  Gainesville. 
Since  January  1999,  he  has  been  with  the  Department  of  Communications 
Engineering,  College  of  Air  Traffic  Management,  Civil  Aviation  University  of 
China,  Tianjin,  China,  where  he  is  currently  a  Chaired  Professor  of  Tianjin. 
His  research  interests  include  spectral  estimation,  feature  extraction  and  image 
formation,  space-time  adaptive  processing,  and  adaptive  arrays  and  their 
applications  to  radar  and  wireless  communications  systems. 


Louis  N.  Cattafesta,  III,  received  the  B.S.  degree 
in  mechanical  engineering  with  Highest  Distinction 
from  the  Pennsylvania  State  University,  University 
Park,  in  1986,  the  M.S.  degree  in  Aeronautics  from 
the  Massachusetts  Institute  of  Technology,  Cam¬ 
bridge,  in  1988,  and  the  Ph.D.  degree  in  mechanical 
engineering  from  the  Pennsylvania  State  University 
in  1992. 

He  is  currently  an  Associate  Professor  with  the 
Department  of  Mechanical  and  Aerospace  Engi¬ 
neering,  University  of  Florida  (UF),  Gainesville. 
Prior  to  joining  UF  in  April  of  1999,  he  was  a  Senior  Research  Scientist  at 
High  Technology  Corporation,  Hampton,  VA,  where  he  was  the  Group  Leader 
of  the  Experimental  and  Instrumentation  Group.  In  1992,  he  joined  High 
Technology  Corporation  as  a  Research  Scientist  at  NASA  Langley  Research 
Center.  His  research  at  NASA  Langley  focused  on  supersonic  laminar  flow 
control  and  pressure-  and  temperature- sensitive  paint  measurement  techniques. 
He  subsequently  became  involved  in  active  control  of  flow-induced  cavity 
oscillations.  His  current  research  interests  lie  in  active  flow  and  noise  control, 
particularly  the  modeling  and  design  of  piezoelectric  actuators  and  MEMS 
sensors  and  also  the  development  and  implementation  of  real-time,  adaptive 
flow  control  schemes. 

Dr.  Cattafesta  is  a  member  of  the  AIAA  Aerodynamic  Measurement  Tech¬ 
nology  Technical  Committee. 


Srihari  Rasmi  received  the  B.Tech.  degree  in  naval 
architecture  from  the  Indian  Institute  of  Technology, 
Madras,  India  in  1999.  He  is  currently  pursuing 
the  M.S  degree  in  the  areas  of  experimental  fluid 
dynamics  and  signal  processing  at  the  Department 
of  Mechanical  and  Aerospace  Engineering  at  the 
University  of  Florida,  Gainesville. 

From  December  1998  to  April  1999,  he  received 
a  scholarship  to  pursue  his  research  interests  related 
to  stability  of  ships  at  Hochshule  Bremen  in  Bremen, 
Germany.  From  August  to  December  of  1999,  he  was 
an  Academic  Visitor  at  the  Ship  Stability  Research  Centre,  University  of  Strath¬ 
clyde,  Glasgow,  U.K. 


Mark  Sheplak  received  the  B.S.,  M.S.,  and  Ph.D.  de¬ 
grees  in  mechanical  engineering  from  Syracuse  Uni¬ 
versity,  Syracuse,  NY,  in  1989,  1992,  and  1995,  re¬ 
spectively. 

He  is  currently  an  Associate  Professor  with 
the  Department  of  Mechanical  and  Aerospace 
Engineering,  University  of  Florida  (UF).  Prior  to 
joining  UF  in  1998,  he  was  a  Postdoctoral  Associate 
at  the  Massachusetts  Institute  of  Technology’s 
Microsystems  Technology  Laboratories,  Cambridge, 
from  1995  to  1998.  During  his  Ph.D.  studies  he 
was  a  GSRP  Fellow  at  NASA-LaRC,  Hampton,  VA,  from  1992  to  1995.  His 
current  research  focuses  on  the  design,  fabrication,  and  characterization  of 
high-performance,  instrumentation-grade,  MEMS -based  sensors  and  actuators 
that  enable  the  measurement,  modeling,  and  control  of  various  physical 
properties. 


